C*************************************************************
C*                    PROGRAM : GRID3O                       *
C*-----------------------------------------------------------*
C*DATE WHEN LAST MODIFICATIONS INTRODUCED: SEPTEMBER 29,1981.*
C*-----------------------------------------------------------*
C*   COMPUTER PROGRAM FOR FAST GENERATION OF 3-D BOUNDARY    *
C*     CONFORMING PERIODIC MULTILEVEL "O"-TYPE GRIDS FOR     *
C*   1. AXIAL TURBOMACHINERY STATOR                          *
C*   2. AXIAL TURBOMACHINERY ROTOR                           *
C*   3. HELICOPTER ROTOR IN HOVER                            *
C*   4. PROPELLER-TYPE WIND TURBINE ROTOR                    *
C*   5. PROPELLER OR PROP-FAN IN FREE AIR                    *
C*   6. PROPELLER OR PROP-FAN IN AXISYMMETRIC TUNNEL         *
C*   7. FINNED MISSILE OR WING-BODY IN FREE AIR              *
C*   8. FINNED MISSILE OR WING-BODY IN AXISYMMETRIC TUNNEL   *
C*   9. MARINE PROPELLER                                     *
C*                           -*-                             *
C*                DR. DJORDJE S. DULIKRAVICH                 *
C*************************************************************
C*    NASA LEWIS RESEARCH CENTER, CLEVELAND, OHIO 44135      *
C* COMPUT.FLUID MECH.BRANCH,5-9,PHONE:(216)433-4000,EXT.6859 *
C*************************************************************
C---   KX    = NO. OF GRID CELLS ON THE AIRFOIL (ON COARSE GRID)
C---   KY    = NO. OF ELLIPTIC LAYERS AROUND BLADE(ON COARSE GRID)
C---   KZ    = NO. OF GRID CELLS ALONG THE BLADE(ON COARSE GRID)
C---   KT    = NO. OF GRID CELLS FROM HUB TO TIP(ON COARSE GRID)
C---   NM    = NO. OF GRIDS TO BE GENERATED
C---           NM>0 : (X,Y,Z) COORDINATES IN OUTPUT
C---           NM<0 : (X,THETA,R) COORDINATES IN OUTPUT
C---   ND    = NO. OF INPUT POINTS ON THE DUCT
C---   NH    = NO. OF INPUT POINTS ON THE HUB (NH .GE. ND)
C---   NB    = NO. OF BLADES
C---   SA    = BLADE SETTING ANGLE (DEGREES)
C---   DZ    = LENGTH UNIT CONVERSION FACTOR
C---   RT    = ROTOR RADIUS (OR WING HALF-SPAN), (LENGTH).
C---   XHUB  = X DISTANCE OF INPUT PNTS ON HUB, (LENGTH).
C---   XDUCT = X DISTANCE OF INPUT PNTS ON DUCT, (LENGTH).
C---   RHUB  = RADIUS OF THE HUB, (LENGTH).
C---   RDUCT = RADIUS OF THE DUCT, (LENGTH).
C---   XLEAD = X COORD. OF THE BLADE L.E. ON N-TH INPUT PLANE, (LENGTH).
C---   YLEAD = Y COORD. OF THE BLADE L.E. ON N-TH INPUT PLANE, (LENGTH).
C---   ZLEAD = Z COORD. OF THE N-TH INPUT PLANE, (LENGTH).
C---   CHORD = MAX. BLADE CHORD LENGTH ON N-TH INPUT PLANE, (LENGTH).
C---   RO1   = RADIUS OF THE L.E./CHORD ON N-TH INPUT PLANE
C---   RO2   = RADIUS OF THE T.E./CHORD ON N-TH INPUT PLANE
C---   TWIST = BLADE TWIST ANGLE ON THE N-TH INPUT PLANE (DEGREES)
C---           TWIST>0 : NOSE DOWN
C---           TWIST<0 : NOSE UP
C---   MAXP  = NUMBER OF INPUT POINTS ON N-TH INPUT PLANE
C --------------------------------------------------------------------
C---NOTE: MAX. NO. OF GRID CELLS WITH PRESENT COMMON REGION:KX=192;KY=48
C---NOTE: MAX. NO. OF INPUT PNTS. AT ANY STATION IS: MAXP=199
C---NOTE: MAX. NO. OF INPUT PLANES IS (51)
C---NOTE: MAX. NO. OF INPUT POINTS ON HUB OR SHROUD IS (51)
C ---------------------------------------------------------------------
C--- NOTE: PROGRAM "GRID3O" CONSISTS OF THIS "MAIN" PROGRAM AND
C---      THE FOLLOWING SUBROUTINES: CENTRL,CONMAP,REMAP,XYZINF,
C---      SPLIF,INTPL,SURF,PISMO.
C-----------------------------------------------------------------------
      COMMON/BLK1/ HR(200),VR(200),XHUB(51),RHUB(51),XDUCT(51)
     1            ,RDUCT(51),XCELL,YCELL,ZCELL,PI,PITCH,MAXP,IMID,ILHS
     2            ,IRHS,MAXXM,MAXX,MAXXP,MAXYM,MAXY,MAXYP,MAXZM,MAXZ
     3            ,MAXZP,KTIP,MHUB,MDUCT
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      COMMON/BLK3/ X(200,51),Y(200,51),Z(200,51)
      COMMON/BLK4/ XP(200,51),YP(200,51),ZP(200,51)
      DOUBLE PRECISION XCELL,YCELL,ZCELL,REST,BLADE,PITCH,TH,RHE,
     1DABS,PI,RADIAN,XDUCT,RDUCT,ZLEAD,YLEAD,TWIST,XLEAD,RO1,RO2,
     2HR,VR,PP,QQ,DSIN,DCOS,TETAS,S,F,SI,FI,FP,FPP,FPPP,YOR,XOR,DVAPI,
     3DSQRT,DATAN,DATAN2,RO,R,RSF,XHUB,RHUB,SL1,SL2,RB,BCELL,XLE,YLE,
     4DTAN,SA,DARSIN,RAD,RH,XP,YP,ZOR,RT,RBTIP,CALE,AA,BB,DXD,R2,ANGL
      WRITE(6,1)
1     FORMAT(1H1,///,27X,43H3-D MULTILEVEL BOUNDARY CONFORMING PERIODIC,
     1  /,38X,24H'O'-TYPE GRID GENERATION)
      WRITE(6,2)
2     FORMAT(26X,46HN.A.S.A. LEWIS RESEARCH CENTER - C.F.M. BRANCH)
      WRITE(6,3)
3     FORMAT(26X,46H(GRID3O) DEVELOPED BY : DJORDJE S. DULIKRAVICH,/,
     1   26X,46(1H=))
      WRITE(6,4)
4     FORMAT(/,43X,16HINPUT PARAMETERS)
5     FORMAT(20A4)
6     FORMAT(7(3X,I2),3X,I3,3X,F6.2,2(3X,F8.4))
C------ READING INPUT DATA
      READ(5,5) (RAD(I),I=1,20)
      READ(5,6) KX,KY,KZ,KT,NM,NH,ND,NB,SA,DZ,RT
      WRITE(6,7) (RAD(I),I=1,20)
7     FORMAT(/,25X,20A4)
      WRITE(6,9)
9     FORMAT(18X,2HKX,5X,2HKY,5X,2HKZ,5X,2HKT,5X,2HNM,5X,2HNB,5X,
     12HSA,5X,2HDZ,7X,2HRT)
      WRITE(6,10) KX,KY,KZ,KT,NM,NB,SA,DZ,RT
10    FORMAT(18X,I2,5X,I2,5X,I2,5X,I2,5X,I2,4X,I3,3X,F6.2,2(2X,F7.4))
      IF(KZ .EQ. KT) GO TO 15
      IF((KZ-KT) .GE. 3) GO TO 15
      WRITE(6,11)
11    FORMAT(///,25X,34HERROR : INCREASE 'KZ' IN THE INPUT,//)
      STOP
15    CONTINUE
      RADIAN=57.2957795130823D+00
      PI=3.14159265358979D+00
      DVAPI=2.0D+00*PI
      CALE=0.4D+00
      LXYZ=0
      IF(NM .GT. 0) LXYZ=1
      NM=IABS(NM)
      BLADE=NB*1.0D+00
      PITCH=DVAPI/BLADE
      NAM=1
      IF(NM .EQ. 2) NAM=2
      IF(NM .EQ. 3) NAM=4
      KX=KX*NAM
      KY=KY*NAM
      KZ=KZ*NAM
      KT=KT*NAM
      XCELL=KX*1.0D+00
      YCELL=KY*1.0D+00
      ZCELL=KZ*1.0D+00
      BCELL=KT*1.0D+00
      MAXX=KX+2
      MAXXM=MAXX-1
      MAXXP=MAXX+1
      MAXY=KY+2
      MAXYM=MAXY-1
      MAXYP=MAXY+1
      MAXZ=KZ+2
      MAXZM=MAXZ-1
      MAXZP=MAXZ+1
      KTIP=KT+2
      KTIPM=KTIP-1
      KTIPP=KTIP+1
      IMID=(MAXX+2)/2
      ILHS=IMID-1
      IRHS=IMID+1
      MHUB=NH
      MDUCT=ND
      MRAD=MAX0(MHUB,MDUCT)
C------ READING OF XHUB,RHUB,XDUCT,RDUCT
18    FORMAT(4F10.6)
      READ(5,18) (XHUB(I),RHUB(I),XDUCT(I),RDUCT(I),I=1,MRAD)
      WRITE(6,21)
21    FORMAT(/,17X,7HXHUB(I),10X,7HRHUB(I),10X,8HXDUCT(I),
     110X,8HRDUCT(I),10X,1HI)
      WRITE(6,22) (XHUB(I),RHUB(I),XDUCT(I),RDUCT(I),I,I=1,MRAD)
22    FORMAT(10X,4E17.6,I10)
      WRITE(6,24)
24    FORMAT(/,10X,5HXLEAD,9X,5HZLEAD,9X,5HYLEAD,9X,5HCHORD,
     18X,3HRO1,5X,3HRO2,3X,5HTWIST,3X,1HN)
C------ THE LARGE INPUT READING LOOP STARTS HERE
C------ NF=1 MEANS THAT THE AIRFOIL SHAPE AT THIS INPUT PLANE IS
C------ DIFFERENT FROM THE SHAPE ON THE PREVIOUS INPUT PLANE
C------ NF=0 MEANS THAT AIRFOIL SHAPE AT THIS INPUT PLANE IS
C------ THE SAME AS AT THE  PREVIOUS INPUT PLANE
C------ NF=2 MEANS THAT THIS IS THE LAST INPUT PLANE AND THAT THE
C------ LOCAL AIRFOIL SHAPE IS SAME AS AT THE PREVIOUS INPUT PLANE
C------ NF=3 MEANS THAT THIS IS THE LAST INPUT PLANE AND THAT THE
C------ LOCAL SHAPE IS NOT THE SAME AS AT THE PREVIOUS INPUT PLANE
      NP=0
      NF=1
      N=1
25    CONTINUE
      READ(5,28) XLEAD,ZLEAD,YLEAD,CHORD,RO1,RO2,TWIST,NF,MAXP
28    FORMAT(7F10.6,2I3)
      IF(ZLEAD .LE. RT) GTC=PITCH*RT/CHORD
      IF(ZLEAD .GT. RT) CHORD=PITCH*ZLEAD/GTC
      WRITE(6,30) XLEAD,ZLEAD,YLEAD,CHORD,RO1,RO2,TWIST,N
30    FORMAT(5X,4E14.6,2F8.4,F7.3,I4)
      IF(NF .EQ. 0 .OR. NF .EQ. 2) GO TO 100
      M8=(MAXP+7)/8
      P1=MAXP*1.00001
      I1=P1/2
      P2=2*I1+0.2
      ISYMM=0
      IF(P1 .GT. P2) ISYMM=1
      IF(ISYMM .EQ. 0) MAXP=2*MAXP-1
C------ READ (X,Y) INPUT PAIRS OF AIRFOIL COORDINATES
35    FORMAT(8F10.7)
      I1=1
      I2=8
      DO 40 II=1,M8
      READ(5,35) (HR(I),I=I1,I2)
      I1=I2+1
      I2=I1+7
40    CONTINUE
      I1=1
      I2=8
      DO 45 II=1,M8
      READ(5,35) (VR(I),I=I1,I2)
      I1=I2+1
      I2=I1+7
45    CONTINUE
      MPHF=(MAXP+1)/2
      IF(ISYMM .EQ. 1) GO TO 50
C------ THE CASE OF A SYMMETRIC AIRFOIL
      DO 47 II=1,MPHF
      I=MAXP+1-II
      HR(I)=HR(II)
      VR(I)=-VR(II)
47    CONTINUE
50    CONTINUE
C------ INPUT BLADE CHORD AND L.E. POINT
      HR(1)=0.5D+00*(HR(1)+HR(MAXP))
      VR(1)=0.5D+00*(VR(1)+VR(MAXP))
      HR(MAXP)=HR(1)
      VR(MAXP)=VR(1)
      CH=-1.0D+00
      DO 55 I=8,MAXP
      ZOR=CH
      XOR=HR(I)-HR(1)
      YOR=VR(I)-VR(1)
      CH=XOR*XOR+YOR*YOR
      IF(CH .LT. ZOR) GO TO 60
55    CONTINUE
60    CONTINUE
      ILE=I-1
      CH=DSQRT(ZOR)
C------ BLADE COORDINATE NORMALIZATION AND PRE-ROTATION
      DO 65 I=1,MAXP
      HR(I)=HR(I)/CH
      VR(I)=VR(I)/CH
65    CONTINUE
      XLE=HR(ILE)
      YLE=VR(ILE)
      PP=VR(1)-YLE
      QQ=HR(1)-XLE
      DO 70 I=1,MAXP
      XOR=HR(I)-XLE
      YOR=VR(I)-YLE
      HR(I)=QQ*XOR+PP*YOR
      VR(I)=QQ*YOR-PP*XOR
70    CONTINUE
C------ EQUAL NO. OF INPUT POINTS ON ALL INPUT PLANES
      S(1)=0.0D+00
      F(1)=0.0D+00
      DO 72 I=2,MAXP
      XOR=HR(I)-HR(I-1)
      YOR=VR(I)-VR(I-1)
      R2=XOR*XOR+YOR*YOR
      S(I)=S(I-1)+DSQRT(R2)
      II=MAXP+2-I
      XOR=HR(II)-HR(II-1)
      YOR=VR(II)-VR(II-1)
      R2=XOR*XOR+YOR*YOR
      F(I)=F(I-1)+DSQRT(R2)
72    CONTINUE
      AB=0.5D+00*(S(MAXP)+F(MAXP))
      AA=0.5D+00*(S(ILE)+F(MAXP)-F(ILE))
      BB=AB-AA
      A=AA/S(ILE)
      B=S(ILE)
      DO 75 I=1,ILE
      S(I)=S(I)*A
75    CONTINUE
      A=B
      B=BB/(S(MAXP)-A)
      ILEP=ILE+1
      DO 78 I=ILEP,MAXP
      S(I)=AA+B*(S(I)-A)
78    CONTINUE
      DXD=(IMID-2)*1.0D+00
      DO 80 I=2,IMID
      II=I-2
      DX=II/DXD
      XOR=DVAPI*DX
      DX=DX-0.10D+00*DSIN(XOR)
      SI(I)=AA*DX
      SI(IMID+II)=AA+BB*DX
80    CONTINUE
C------ CUBIC SPLINE FITTING AND INTERPOLATION
      DO 82 I=1,MAXP
      F(I)=HR(I)
82    CONTINUE
      CALL SPLIF(1,MAXP)
      CALL INTPL(2,MAXX,1,MAXP)
      DO 85 I=2,MAXX
      HR(I)=FI(I)
85    CONTINUE
      DO 90 I=1,MAXP
      F(I)=VR(I)
90    CONTINUE
      CALL SPLIF(1,MAXP)
      CALL INTPL(2,MAXX,1,MAXP)
      DO 95 I=2,MAXX
      VR(I)=FI(I)
95    CONTINUE
      HR(2)=0.5D+00*(HR(2)+HR(MAXX))
      VR(2)=0.5D+00*(VR(2)+VR(MAXX))
      HR(MAXX)=HR(2)
      VR(MAXX)=VR(2)
100   CONTINUE
      TWIST=(TWIST+SA)/RADIAN
      PP=DSIN(TWIST)
      QQ=DCOS(TWIST)
      RAD1(N)=RO1
      RAD2(N)=RO2
      RO=ZLEAD*ZLEAD
C------ ROTATION OF THE AIRFOIL TO THE CORRECT INCLINATION
      DO 120 I=2,MAXX
      SI(I)=(HR(I)*QQ-VR(I)*PP)*CHORD+XLEAD
      YOR=(VR(I)*QQ+HR(I)*PP)*CHORD+YLEAD
      R2=YOR*YOR+RO
      RAD(I)=DSQRT(R2)
      YOR=YOR/ZLEAD
      YP(I,N)=DATAN(YOR)
120   CONTINUE
C------ SI(I) : AXIAL COORD.
C------ RAD(I): RADIAL COORD.
C------ YP(I) : ANGULAR COORD.
C------ THE MOST RIGHTWARD PNT ON THE AIRFOIL SURFACE
      MXM=MAXX-5
      DO 130 I=MXM,MAXX
      II=MAXX+2-I
      IF(SI(I) .LE. SI(I-1)) GO TO 140
      IF(SI(II) .LE. SI(II+1)) GO TO 150
130   CONTINUE
      ITRAIL=MAXX
      GO TO 160
140   CONTINUE
      ITRAIL=I-1
      GO TO 160
150   CONTINUE
      ITRAIL=II+1
160   CONTINUE
C------ THE MOST LEFTWARD PNT ON THE AIRFOIL SURFACE
      DO 170 I=5,MXM
      IF(SI(I) .GE. SI(I-1)) GO TO 180
170   CONTINUE
180   CONTINUE
      ILEAD=I-1
C------ INTERPOLATION OF (RHUB) AT X-AIRFOIL
      DO 190 I=1,MHUB
      S(I)=XHUB(I)
      F(I)=RHUB(I)
190   CONTINUE
      CALL SPLIF(1,MHUB)
      IF(ITRAIL .GE. 10) GO TO 200
      CALL INTPL(ITRAIL,ILEAD,1,MHUB)
      CALL INTPL(ILEAD,MAXX,1,MHUB)
      IF(ITRAIL .GT. 1) CALL INTPL(2,ITRAIL,1,MHUB)
      GO TO 220
200   CONTINUE
      CALL INTPL(2,ILEAD,1,MHUB)
      CALL INTPL(ILEAD,ITRAIL,1,MHUB)
      IF(MAXX .GT. ITRAIL) CALL INTPL(ITRAIL,MAXX,1,MHUB)
220   CONTINUE
      DO 240 I=2,MAXX
      RH(I)=FI(I)
240   CONTINUE
C------ INTERPOLATION OF (RDUCT) AT X-AIRFOIL
      DO 260 I=1,MDUCT
      S(I)=XDUCT(I)
      F(I)=RDUCT(I)
260   CONTINUE
      CALL SPLIF(1,MDUCT)
      IF(ITRAIL .GE. 10) GO TO 280
      CALL INTPL(ITRAIL,ILEAD,1,MDUCT)
      CALL INTPL(ILEAD,MAXX,1,MDUCT)
      IF(ITRAIL .GT. 1) CALL INTPL(2,ITRAIL,1,MDUCT)
      GO TO 300
280   CONTINUE
      CALL INTPL(2,ILEAD,1,MDUCT)
      CALL INTPL(ILEAD,ITRAIL,1,MDUCT)
      IF(MAXX .GT. ITRAIL) CALL INTPL(ITRAIL,MAXX,1,MDUCT)
300   CONTINUE
C------ INTERPOLATED (XBAR) AND (RBAR) VALUES
      DO 320 I=2,MAXX
      XP(I,N)=SI(I)
      RO=RAD(I)-RH(I)
      ZP(I,N)=CALE+RO/(FI(I)-RH(I))
320   CONTINUE
      NP=N
      N=N+1
      IF(NF .LT. 2) GO TO 25
340   CONTINUE
C------ RADIAL COORDINATE SHEARING
C------ NEW RADIAL COORD. : HUB=CALE ; DUCT=CALE+1. ; TIP=RBTIP
C------ MAXIMUM VALUE FOR (SL1) IS 0.15
      SL1=0.06D+00
      REST=(RT-RH(2))/(FI(2)-RH(2))
      DO 360 K=2,KTIP
      R=((K-2)*1.0D+00)/BCELL
      ZOR=DVAPI*R
      SI(K)=CALE+REST*(R-SL1*DSIN(ZOR))
360   CONTINUE
      IF(MAXZ .EQ. KTIP) GO TO 380
      SL1=REST*(1.0D+00+DVAPI*SL1)
      REST=1.0D+00-REST
      SL1=(1.0D+00-SL1/REST)/PI
      RO=ZCELL-BCELL
      DO 370 K=KTIP,MAXZ
      R=((K-KTIP)*1.0D+00)/RO
      ZOR=PI*R
      SI(K)=SI(KTIP)+REST*(R-SL1*DSIN(ZOR))
370   CONTINUE
380   CONTINUE
      DO 480 K=2,MAXZ
      RAD(K)=SI(K)
480   CONTINUE
C------ LEADING EDGE RADII
      DO 500 N=1,NP
      S(N)=ZP(ILEAD,N)
      F(N)=RAD1(N)
500   CONTINUE
      CALL SPLIF(1,NP)
      CALL INTPL(2,MAXZ,1,NP)
      DO 520 K=2,MAXZ
      RAD1(K)=FI(K)
520   CONTINUE
C------ TRAILING EDGE RADII
      DO 540 N=1,NP
      S(N)=ZP(ITRAIL,N)
      F(N)=RAD2(N)
540   CONTINUE
      CALL SPLIF(1,NP)
      CALL INTPL(2,MAXZ,1,NP)
      DO 560 K=2,MAXZ
      RAD2(K)=FI(K)
560   CONTINUE
      DO 675 I=2,MAXX
C------ SPANWISE INTERPOLATION OF X-COORDINATE
      DO 580 N=1,NP
      S(N)=ZP(I,N)
      F(N)=XP(I,N)
580   CONTINUE
      CALL SPLIF(1,NP)
      CALL INTPL(2,MAXZ,1,NP)
      DO 600 K=2,MAXZ
      XP(I,K)=FI(K)
600   CONTINUE
C------ SPANWISE INTERPOLATION OF (THETA)-COORD.
      DO 625 N=1,NP
      F(N)=YP(I,N)
625   CONTINUE
      CALL SPLIF(1,NP)
      CALL INTPL(2,MAXZ,1,NP)
      DO 650 K=2,MAXZ
      YP(I,K)=FI(K)
650   CONTINUE
675   CONTINUE
C------ NOW WE KNOW X,THETA AND RBAR (WHICH IS CONST. AT
C------ EACH K-TH POSITION). THESE SHAPES ARE ON THE CYLIND.
C------ RBAR SURFACES I.E., ON THE AXISYMMETRIC RBAR=R(X) SURF.
C------ GEOMETRY WILL NOW BE TREATED AS IF THE HUB AND THE
C------ DUCT ARE TWO COAXIAL CIRCULAR CYLINDERS
      WRITE(6,700)
700   FORMAT(//,40X,17HOUTPUT PARAMETERS)
      WRITE(6,725)
725   FORMAT(/,10X,5HXLEAD,9X,5HZLEAD,9X,5HYLEAD,9X,5HCHORD,
     18X,3HRO1,5X,3HRO2,3X,5HTWIST,3X,1HK)
      DO 950 K=2,MAXZ
      N1=NM+10
      RBAR=RAD(K)
      DO 750 I=2,MAXX
      S(I)=XP(I,K)
      F(I)=YP(I,K)*RBAR
750   CONTINUE
      S(2)=0.5D+00*(S(2)+S(MAXX))
      F(2)=0.5D+00*(F(2)+F(MAXX))
      F(MAXX)=F(2)
      S(MAXX)=S(2)
      CALL CENTRL(K)
      CALL SURF(K,CALE)
      CALL XYZINF(K)
C------ LENGTH UNITS CONVERSION
      DO 775 I=1,MAXXP
      DO 775 J=2,MAXY
      X(I,J)=X(I,J)*DZ
      Z(I,J)=Z(I,J)*DZ
775   CONTINUE
      CALL PISMO(K,N1,LXYZ)
      DO 800 I=5,MAXX
      IF(X(I,MAXY) .GT. X(I-1,MAXY)) GO TO 825
800   CONTINUE
825   IL=I-1
      XL=X(IL,MAXY)
      RL=Z(IL,MAXY)
      TL=Y(IL,MAXY)
      ZL=RL*COS(TL)
      YL=RL*SIN(TL)
      XT=X(2,MAXY)
      RT=Z(2,MAXY)
      TT=Y(2,MAXY)
      ZT=RT*COS(TT)
      YT=RT*SIN(TT)
      CH=(XT-XL)*(XT-XL)+(YT-YL)*(YT-YL)
      CH=SQRT(CH)
      TW=(YT-YL)/CH
      TW=180.*ARSIN(TW)/PI
      WRITE(6,30) XL,ZL,YL,CH,RAD1(K),RAD2(K),TW,K
      LOOP=1
850   CONTINUE
      N1=N1-1
      IF(N1 .EQ. 10) GO TO 925
      K1=K/2
      K2=K1*2
      IF(K2 .EQ. K) GO TO 875
      GO TO 950
875   CONTINUE
      K3=(K+2)/4
      K4=K3*4
      IF(K4 .EQ. K .AND. LOOP .EQ. 2) GO TO 925
      II=1
      DO 900 I=2,MAXX,2
      II=II+1
      JJ=1
      DO 900 J=2,MAXY,2
      JJ=JJ+1
      X(II,JJ)=X(I,J)
      Y(II,JJ)=Y(I,J)
      Z(II,JJ)=Z(I,J)
900   CONTINUE
      MAXX=(MAXX+2)/2
      MAXY=(MAXY+2)/2
      MAXXP=MAXX+1
      MAXYP=MAXY+1
      MAXXM=MAXX-1
      MAXYM=MAXY-1
      IMID=(MAXX+2)/2
      IRHS=IMID+1
      ILHS=IMID-1
      CALL XYZINF(K)
      CALL PISMO(K,N1,LXYZ)
      LOOP=LOOP+1
      GO TO 850
925   CONTINUE
      MAXX=KX+2
      MAXY=KY+2
      MAXXM=MAXX-1
      MAXXP=MAXX+1
      MAXYM=MAXY-1
      MAXYP=MAXY+1
      IMID=(MAXX+2)/2
      IRHS=IMID+1
      ILHS=IMID-1
950   CONTINUE
      MX=KX+3
      MY=KY+3
      MZ=KZ+1
      MT=10+NM
      IF(LXYZ .EQ. 1) WRITE(6,975) MX,MY,MZ,MT
      IF(LXYZ .EQ. 0) WRITE(6,330) MX,MY,MZ,MT
975   FORMAT(//,15X,24H(X,Y,Z) COORDINATES OF (,I3,1H*,I2,1H*,I2,1H),
     1   1X,23HPOINTS WRITTEN ON UNIT ,I2,/)
330   FORMAT(//,15X,24H(X,T,R) COORDINATES OF (,I3,1H*,I2,1H*,I2,1H),
     1   1X,23HPOINTS WRITTEN ON UNIT ,I2,/)
      IF(MT .EQ.11) GO TO 1000
      MX=KX/2+3
      MY=KY/2+3
      MZ=KZ/2+1
      MT=MT-1
      IF(LXYZ .EQ. 1) WRITE(6,975) MX,MY,MZ,MT
      IF(LXYZ .EQ. 0) WRITE(6,330) MX,MY,MZ,MT
      IF(MT .EQ. 11) GO TO 1000
      MX=KX/4+3
      MY=KY/4+3
      MZ=KZ/4+1
      MT=MT-1
      WRITE(6,975) MX,MY,MZ,MT
1000  CONTINUE
      RETURN
      END
      SUBROUTINE CENTRL(K)
C********************************************************
C* ROUTINE THAT ITERATIVELY DETERMINES SLIT HALF-LENGTH *
C* IN CIRCLE PLANE AND SETS UP GEOMETRY FOR MAPPING AND *
C*        REMAPPING FROM THE COMPUTATIONAL PLANE        *
C********************************************************
      COMMON/BLK1/ HR(200),VR(200),XHUB(51),RHUB(51),XDUCT(51)
     1            ,RDUCT(51),XCELL,YCELL,ZCELL,PI,PITCH,MAXP,IMID,ILHS
     2            ,IRHS,MAXXM,MAXX,MAXXP,MAXYM,MAXY,MAXYP,MAXZM,MAXZ
     3            ,MAXZP,KTIP,MHUB,MDUCT
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      COMMON/BLK3/ X(200,51),Y(200,51),Z(200,51)
      COMMON/BLK4/ XP(200,51),YP(200,51),ZP(200,51)
      DOUBLE PRECISION HR,VR,RHUB,XDUCT,RDUCT,XCELL,YCELL,ZCELL,
     1 PI,PITCH,S,F,SI,FI,FP,FPP,FPPP,RAD,RH,DVAPI,XS,YS,XOR,YOR,
     2 CH,DSQRT,TWIST,AC,BC,DATAN,PP,QQ,DSIN,DCOS,SL,C1C2,ALM,
     3 R1,R2,EMP,EMM,PP2,EM,DEXP,EM2,FO,A,B,DRD,AA,RBAR,XHUB,DABS,
     4 DVAQQ,DVAPP,HALF,XSI,DLOG,DARSIN,EMOLD
      DVAPI=2.0D+00*PI
C------ BLADE CHORD AND L.E. POINT ON (X,RBAR*THETA) PLANE
      CH=-1.0D+00
      DO 5 I=8,MAXX
      ZOR=CH
      XOR=S(I)-S(2)
      YOR=F(I)-F(2)
      CH=XOR*XOR+YOR*YOR
      IF(CH .LT. ZOR) GO TO 6
5     CONTINUE
6     CONTINUE
      MID=I-1
      CH=DSQRT(ZOR)
      RBAR=RAD(K)
      XS=S(MID)
      YS=F(MID)
      XOR=S(2)-XS
      YOR=F(2)-YS
      AC=-YOR/CH
      BC=XOR/CH
      NTW=1
      IF(YOR .GT. 1.0D-05) NTW=-1
      NU=1
      IF(XOR .LT. 0.0D+00) NU=-1
      QQ=BC*NTW
      PP=AC*NTW
C------ STRETCHING FACTOR
      SL=DVAPI*CH/(RBAR*PITCH)
      R1=RAD1(K)/DSQRT(RBAR)
      R2=RAD2(K)/DSQRT(RBAR)
      C1C2=(1.0D+00-(R1+R2)*0.5D+00)*SL
C------ SLIT HALF-LENGTH IN CIRCLE PLANE
      EM=0.25D+00*C1C2*DABS(QQ)
      EM=DEXP(EM)
      EM=(EM-1.0D+00)/(EM+1.0D+00)
      PP2=PP*PP
      DVAQQ=2.0D+00*QQ
      DVAPP=2.0D+00*PP
      DO 20 N=1,100
      EM2=EM*EM
      EMP=1.0D+00+EM2
      EMM=1.0D+00-EM2
      HALF=EM*DVAQQ/EMM
      XSI=1.0D+00+HALF*HALF
      XOR=HALF+DSQRT(XSI)
      XOR=QQ*DLOG(XOR)
      YOR=EM*DVAPP/EMP
      YOR=PP*DARSIN(YOR)
      A=-2.0D+00*(XOR+YOR)+0.5D+00*C1C2
      XOR=EMP/DSQRT(XSI)
      XOR=XOR/(EMM*EMM)
      YOR=EMP*EMP-4.0D+00*EM2*PP2
      YOR=EMM/(EMP*DSQRT(YOR))
      B=-4.0D+00*(QQ*QQ*XOR+PP2*YOR)
      EMOLD=EM
      EM=EMOLD-A/B
      XOR=(EM-EMOLD)/EM
      IF(DABS(XOR) .LT. 1.0D-12) GO TO 30
20    CONTINUE
      WRITE(6,25) EM,XOR
25    FORMAT(//,10X,18HDIVERGES IN CENTRL,5X,3HEM=,E12.5,4HXOR=,E12.5)
      STOP
30    CONTINUE
C------ ANGULAR POSITION OF TRAILING EDGE IN CIRCLE PLANE
      EM2=EM*EM
      FO=(1.0D+00-EM2)/(1.0D+00+EM2)
      FO=FO*PP/QQ
      FO=DATAN(FO)+((NTW+1)/2)*PI
C------ ROTATION OF AIRFOIL INTO HORIZONTAL POSITION
      A=(NU*1.0D+00)/CH
      B=PI*QQ
      ALM=0.5D+00*C1C2+R1*0.5D+00*SL+PI*PP
      DO 35 I=2,MAXX
      XOR=((S(I)-XS)*BC-(F(I)-YS)*AC)*A
      YOR=((F(I)-YS)*BC+(S(I)-XS)*AC)*A
      HR(I)=XOR*SL-ALM
      VR(I)=YOR*SL+B
35    CONTINUE
      VR(MAXX)=0.5D+00*(VR(2)+VR(MAXX))
      VR(2)=VR(MAXX)-1.0D-06
      VR(MAXX)=VR(MAXX)+1.0D-06
      HR(MAXX)=0.5D+00*(HR(2)+HR(MAXX))
      HR(2)=HR(MAXX)
      CALL CONMAP(EM,FO,PP,QQ,NTW)
      XOR=PI-0.5D+00*(S(MAXX)-S(2))
      S(2)=S(2)-XOR*NTW
      S(MAXX)=S(MAXX)+XOR*NTW
      DRD=S(MAXX)-DVAPI
      AA=0.175D+00/RBAR
      DO 40 I=2,MAXX
      DX=((I-2)*DVAPI)/XCELL
      YOR=2.0D+00*DX
      HR(I)=-PI+DX-AA*DSIN(YOR)
40    CONTINUE
      CALL REMAP(AC,BC,XS,YS,EM,DRD,NTW,SL,ALM,CH,NU,K)
      RETURN
      END
      SUBROUTINE CONMAP(EM,FO,PP,QQ,NTW)
C*************************************************************
C* THIS ROUTINE UTILIZES NEWTON-RAPHSON ITERATIVE PROCEDURE  *
C* TO FIND THE INVERSE FUNCTION Z=Z(W) OF THE BASIC REMAPPING*
C*                        FUNCTION W=W(Z)                    *
C*************************************************************
      COMMON/BLK1/ HR(200),VR(200),XHUB(51),RHUB(51),XDUCT(51)
     1            ,RDUCT(51),XCELL,YCELL,ZCELL,PI,PITCH,MAXP,IMID,ILHS
     2            ,IRHS,MAXXM,MAXX,MAXXP,MAXYM,MAXY,MAXYP,MAXZM,MAXZ
     3            ,MAXZP,KTIP,MHUB,MDUCT
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      DOUBLE PRECISION HR,VR,PI,PITCH,XHUB,RHUB,XDUCT,RDUCT,EM2,
     1XSI,ETA,XSI2,ETA2,DVAMX,DVAMY,DVAMXY,XCELL,YCELL,ZCELL
     2,X2Y2,X2Y2M,D1X,D2X,D3X,D4X,D1,D2,D3,D4,DATAN2,R1Y2,R1
     3,R2,R3,R4,PP,QQ,R1R2,DSIN,DCOS,DLOG,A1,A2,A3,A4,ZR,DVAPI
      DOUBLE PRECISION TO,QO,RO2,RO,BE,BE2,P12,UE,VE,UE1,UE2,S,
     1EM,P3,P4,Q3,Q4,DSINH,DCOSH,DSIGN,DVAXY,DSQRT,ZDOL,F,
     2XSINEW,ETANEW,DIFFX,DIFFY,DABS,XSIO,ETAO,DATAN,TRI,SI,
     3ZI,UEC,DX,DY,DARCOS,CETR,FO,DM,FI,FP,FPP,FPPP,RAD,RH
      DVAPI=2.0D+00*PI
      DM=2.0D+00*EM
      EM2=EM*EM
      DX=1.0D-07
      DY=1.0D-07
      MID=ILHS
      DO 5 I=8,MAXX
      IF(HR(I) .GT. HR(I-1)) GO TO 6
5     CONTINUE
6     CONTINUE
      MID=I-1
C------ INITIAL GUESS
      XSI=DCOS(FO)*0.75D+00
      ETA=DSIN(FO)*0.75D+00
      DO 90 II=2,MAXX
      I=II
      IF(NTW .LT. 0) I=MAXX+2-II
      XSIO=HR(I)
      ETAO=VR(I)
      DO 10 N=1,200
      XSI2=XSI*XSI
      ETA2=ETA*ETA
      DVAMX=DM*XSI
      DVAMY=-DM*ETA
      DVAXY=-2.0D+00*XSI*ETA
      DVAMXY=DVAXY*EM2
      X2Y2=XSI2+ETA2
      X2Y2M=-XSI2+ETA2
      D1X=EM2-X2Y2
      D2X=1.0D+00-EM2*X2Y2
      D1=DATAN2(DVAMY,D1X)
      IF(D1 .LT. 0.0D+00) D1=D1+DVAPI
      D2=DATAN2(DVAMY,D2X)
      D3X=-(EM2+X2Y2M)
      D4X=-(1.0D+00+EM2*X2Y2M)
      TRI=D3X*D3X+DVAXY*DVAXY
      CETR=D4X*D4X+DVAMXY*DVAMXY
      R3=DM/DSQRT(TRI)
      R4=DM/DSQRT(CETR)
      D3=DATAN2(DVAXY,D3X)
      D4=DATAN2(DVAMXY,D4X)
      P3=DSIN(D3)
      Q3=DCOS(D3)
      P4=DSIN(D4)
      Q4=DCOS(D4)
      R1Y2=DVAMY*DVAMY
      R1=D1X*D1X+R1Y2
      R2=D2X*D2X+R1Y2
      R1=DSQRT(R1)/(EM2+X2Y2+DVAMX)
      R2=DSQRT(R2)/(1.0D+00+EM2*X2Y2+DVAMX)
      R1R2=R1*R2
      R1R2=DLOG(R1R2)
      A1=QQ*R1R2+PP*(D2-D1)-XSIO
      R1R2=R1/R2
      R1R2=DLOG(R1R2)
      A2=PP*R1R2+QQ*(D1+D2)-ETAO
      A3=R3*(QQ*Q3-PP*P3)+R4*(QQ*Q4+PP*P4)
      A4=R3*(PP*Q3+QQ*P3)-R4*(PP*Q4-QQ*P4)
      ZDOL=A3*A3+A4*A4
      DIFFX=-(A1*A3+A2*A4)/ZDOL
      DIFFY=(A1*A4-A2*A3)/ZDOL
      XSI=XSI+0.75D+00*DIFFX
      ETA=ETA+0.75D+00*DIFFY
      DIFFX=DIFFX/XSI
      IF(DABS(ETA) .GT. 1.0D-10) DIFFY=DIFFY/ETA
      IF(DABS(DIFFX) .LT. DX .AND. DABS(DIFFY) .LT. DY) GO TO 30
10    CONTINUE
      WRITE(6,20) I,N,DIFFX,DIFFY
20    FORMAT(//,10X,2HI=,I3,5X,2HN=,I3,5X,6HDIFFX=,E13.6,
     1  4X,6HDIFFY=,E13.6,/,30X,18HDIVERGES IN CONMAP,///)
      STOP
30    CONTINUE
      IF(DABS(XSI) .LT. 1.0D-10) GO TO 50
      BE=1.0D+00+X2Y2/EM2
      BE2=BE*BE-4.0D+00*XSI2/EM2
      P12=0.5D+00*(BE-DSQRT(BE2))
      P12=DSQRT(P12)
      IF(P12 .GT. 1.0D+00) P12=1.0D+00
      P12=-DSIGN(P12,XSI)
      VE=DARCOS(P12)-DVAPI
      IF(ETA .GT. 0.0D+00) VE=-VE-DVAPI
      IF(ETA .LT. 0.0D+00 .AND. I .GT. MID) VE=VE+DVAPI
      VE=VE+PI
      UE1=XSI/(EM*P12)
      GO TO 60
50    CONTINUE
      VE=0.5D+00*PI*DSIGN(1.0D+00,VE)
      UE1=1.0D+00+X2Y2/EM2
      UE1=DSQRT(UE1)
60    CONTINUE
      UE2=UE1*UE1-1.0D+00
      UE=DABS(UE1)+DSQRT(UE2)
      S(I)=VE
      F(I)=DLOG(UE)
      IF(I .EQ. 2 .OR. I .EQ. MAXX .OR. I .EQ. MID) GO TO 70
      GO TO 80
70    CONTINUE
C------ INITIAL GUESS FOR THE NEXT POINT
      RO=0.8*DSQRT(X2Y2)
      D1=0.1*NTW+DATAN2(ETA,XSI)
      XSI=DCOS(D1)*RO
      ETA=RO*DSIN(D1)
80    CONTINUE
90    CONTINUE
      F(2)=0.5D+00*(F(2)+F(MAXX))
      F(MAXX)=F(2)
      RETURN
      END
      SUBROUTINE REMAP(AC,BC,XS,YS,EM,DRD,NTW,SL,ALM,CH,NU,K)
C*************************************************************
C* ROUTINE THAT MAPS UNIFORM MESH IN COMPUTATIONAL SPACE BACK*
C*                      INTO THE PHYSICAL SPACE              *
C*************************************************************
      COMMON/BLK1/ HR(200),VR(200),XHUB(51),RHUB(51),XDUCT(51)
     1            ,RDUCT(51),XCELL,YCELL,ZCELL,PI,PITCH,MAXP,IMID,ILHS
     2            ,IRHS,MAXXM,MAXX,MAXXP,MAXYM,MAXY,MAXYP,MAXZM,MAXZ
     3            ,MAXZP,KTIP,MHUB,MDUCT
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      COMMON/BLK3/ X(200,51),Y(200,51),Z(200,51)
      COMMON/BLK4/ XP(200,51),YP(200,51),ZP(200,51)
      DOUBLE PRECISION HR,VR,RDUCT,XDUCT,PI,PITCH,ALM,DX2,SIS,XSI,DXP,
     1 XCELL,YCELL,ZCELL,UE,DCOSH,VE,DSQRT,DSINH,US,EM,PP,QQ,SOS,DYP,
     2 DLOG,DATAN2,T1,T2,DSIN,DCOS,UIN,DRD,DVAPI,YOY,EM2,DX1,SL,C,DY,
     3 ETA,R1,R2,DVAMX,DVAMY,XSI2,ETA2,RR,DABS,RHUB,A,ANA,R,XHUB,C2,
     4 AC,BC,XS,YS,CH,DD,XX,YY,S,F,SI,FI,FP,FPPP,FPP,RAD,RH,AKA,USA,AK
      DOUBLE PRECISION RBAR,TE,TE1,TE2,CON,X2,CONST,DSIGN,XP,YP
      DVAPI=2.0D+00*PI
      MAXXMM=MAXXM-1
      PP=AC*NTW
      QQ=BC*NTW
      RBAR=RAD(K)
      TE1=(F(3)-F(2))/(S(3)-S(2))
      TE2=(F(MAXX)-F(MAXXM))/(S(MAXX)-S(MAXXM))
      TE=0.5D+00*(TE1+TE2)
      CON=DRD/F(2)
      CONST=3.60D+00*(CON-TE)/PI
      IF(DABS(CONST) .GT. 0.30D+00) CONST=0.30D+00
      TE1=(YP(2,K)-YP(3,K))/(XP(2,K)-XP(3,K))
      TE2=(YP(MAXX,K)-YP(MAXXM,K))/(XP(MAXX,K)-XP(MAXXM,K))
      TE=0.5D+00*(TE1+TE2)
      CONST=DSIGN(CONST,TE)
      EM2=EM*EM
      CH=CH*NU
      DD=2.0D+00*EM
      NN=1-NTW*NU
      N2N=2*NTW*NU
      CHSL=SL/CH
      SOS=ALM/CHSL
      SIS=PI*QQ/CHSL
      DY=DVAPI/XCELL
      AK=0.05D+00+((K-2)*0.20D+00)/ZCELL
      DO 140 J=2,MAXY
      I1=2
      I2=IMID
      MMM=2
      IMNS=-1
      US=((J-2)*1.0D+00)/YCELL
      USA=1.0D+00-US
      AKA=USA*AK
      ANA=US*PI
      ANA=US+(((K-2)*0.08D+00)/ZCELL)*DSIN(ANA)
      DXP=ANA*PI
      DXP=US+CONST*DSIN(DXP)
      UIN=DRD*DXP+0.5D+00*(1+NTW)*PI*(DXP-1.0D+00)
      UIN=(S(2)*US+HR(2)*USA)-UIN
      DO 135 N=1,2
      IF(J .NE. 2) GO TO 125
      I1=I1+1
      I2=I2-1
125   CONTINUE
      DO 130 I=I1,I2
      VE=S(I)*US+HR(I)*USA-UIN
      DYP=(I-MMM)*DY
C     UE=F(I)*US*(1.0D+00+AKA*DSIN(DYP))
      UE=F(I)*ANA
      XSI=EM*DCOSH(UE)*DCOS(VE)
      ETA=EM*DSINH(UE)*DSIN(VE)
      DVAMX=DD*XSI
      DVAMY=-DD*ETA
      X2Y2=XSI*XSI+ETA*ETA
      DX1=EM2-X2Y2
      DX2=1.0D+00-EM2*X2Y2
      T1=DATAN2(DVAMY,DX1)
      IF(T1 .LT. 0.0D+00) T1=DVAPI+T1
      M=I/IMID
      IF(J .EQ. 2) T1=PI*(NN+N2N*FLOAT(M))
      T2=DATAN2(DVAMY,DX2)
      R1=DX1*DX1+DVAMY*DVAMY
      R1=DSQRT(R1)/(EM2+X2Y2+DVAMX)
      R2=DX2*DX2+DVAMY*DVAMY
      R2=DSQRT(R2)/(1.0D+00+EM2*X2Y2+DVAMX)
      RR=R1*R2
      RR=DABS(RR)
      XX=QQ*DLOG(RR)+PP*(T2-T1)
      RR=R1/R2
      RR=DABS(RR)
      YY=PP*DLOG(RR)+QQ*(T2+T1)
      XOR=XX/CHSL+SOS
      YOR=YY/CHSL-SIS
C ----- X(I,J) = X COORDINATE
C------ Y(I,J) = THETA COORDINATE
      X(I,J)=XS+(BC*XOR+AC*YOR)
      Y(I,J)=(YS+BC*YOR-AC*XOR)/RBAR
130   CONTINUE
      I1=IMID
      I2=MAXX
      MMM=IMID
      IMNS=1
135   CONTINUE
140   CONTINUE
      RETURN
      END
      SUBROUTINE XYZINF(K)
C*************************************************************
C*    SUBROUTINE THAT DETERMINES(X,THETA,R) VALUES OF THE    *
C*       IMAGINARY POINTS AND POINTS AT AXIAL INFINITIES     *
C*************************************************************
      COMMON/BLK1/ HR(200),VR(200),XHUB(51),RHUB(51),XDUCT(51)
     1            ,RDUCT(51),XCELL,YCELL,ZCELL,PI,PITCH,MAXP,IMID,ILHS
     2            ,IRHS,MAXXM,MAXX,MAXXP,MAXYM,MAXY,MAXYP,MAXZM,MAXZ
     3            ,MAXZP,KTIP,MHUB,MDUCT
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      COMMON/BLK3/ X(200,51),Y(200,51),Z(200,51)
      DOUBLE PRECISION HR,VR,XHUB,RHUB,XDUCT,RDUCT,PI,PITCH,P2,P4,A,
     1XCELL,YCELL,ZCELL,S,F,SI,FI,FP,FPP,FPPP,RAD,RH,XC,X1,X2,XXP,DSIN
     2,DMIN1,DMAX1
C------ X(I,J) : AXIAL COORDINATE
C------ Y(I,J) : THETA-COORDINATE
C------ Z(I,J) : RADIAL COORDINATE
C------ DELIBERATE MODIFICATIONS - POINTS AT INFINITIES
      P2=0.5D+00*PITCH
      P4=0.25D+00*PITCH
      ILHSM=ILHS-1
      IRHSP=IRHS+1
      MAXXMM=MAXXM-1
      M2=MAXXP+1
      MAXX2=MAXX+2
      IAO=MAXX-5
      DO 10 I=8,IAO
      IF(X(I,3) .GT. X(I-1,3)) GO TO 20
10    CONTINUE
20    CONTINUE
C------ UPSTREAM INFINITY
      ILEFT=I-1
      X1=X(ILEFT,3)-(X(ILEFT,4)-X(ILEFT,3))
      X2=0.5D+00*(X(ILHSM,2)+X(ILHSM-1,2))
      X1=DMIN1(X1,X2)
      X(IMID,2)=X1
      X(ILHS,2)=X1
      X(IRHS,2)=X1
      X(ILHSM,2)=X1
      X(IRHSP,2)=X1
      X(IRHSP+1,2)=0.5D+00*(X(IRHSP,2)+X(IRHSP+2,2))
      X(ILHSM-1,2)=X(IRHSP+1,2)
      A=2.0D+00*Y(ILHSM,2)-Y(ILHSM-1,2)
      B=2.0D+00*Y(IRHSP,2)-Y(IRHSP+1,2)
      Y(IMID,2)=0.5D+00*(A+B)
      Y(ILHS,2)=Y(IMID,2)-P4
      Y(IRHS,2)=Y(IMID,2)+P4
      Y(ILHSM,2)=Y(IMID,2)-P2
      Y(IRHSP,2)=Y(ILHSM,2)+PITCH
      Y(ILHSM-1,2)=0.5D+00*(Y(ILHSM-2,2)+Y(ILHSM,2))
      Y(IRHSP+1,2)=Y(ILHSM-1,2)+PITCH
      Z(IMID,2)=2.0D+00*Z(IMID,3)-Z(IMID,4)
      Z(ILHS,2)=Z(IMID,2)
      Z(IRHS,2)=Z(IMID,2)
      Z(ILHSM,2)=Z(IMID,2)
      Z(IRHSP,2)=Z(IMID,2)
      Z(IRHSP+1,2)=0.5D+00*(Z(IRHSP,2)+Z(IRHSP+2,2))
      Z(ILHSM-1,2)=Z(IRHSP+1,2)
      DO 30 I=IAO,MAXXP
      IM=I-1
      IF(X(I,3) .LT. X(IM,3)) GO TO 40
30    CONTINUE
40    CONTINUE
C------ DOWNSTREAM INFINITY
      IRIGHT=IM
      DO 50 I=1,10
      II=12-I
      IF(X(II,3) .GT. X(IRIGHT,3)) IRIGHT=II
50    CONTINUE
      X2=X(IRIGHT,3)+(X(IRIGHT,3)-X(IRIGHT,4))
      X1=0.5D+00*(X(4,2)+X(5,2))
      X2=DMAX1(X1,X2)
      X(2,2)=X2
      X(MAXX,2)=X2
      X(3,2)=X2
      X(MAXXM,2)=X2
      X(4,2)=X2
      X(MAXXMM,2)=X2
      X(5,2)=0.5D+00*(X(4,2)+X(6,2))
      X(MAXXMM-1,2)=X(5,2)
      Y(2,2)=Y(2,3)
      Y(MAXX,2)=Y(2,2)
      Y(3,2)=Y(2,2)-P4
      Y(MAXXM,2)=Y(2,2)+P4
      Y(4,2)=Y(2,2)-P2
      Y(MAXXMM,2)=Y(4,2)+PITCH
      Y(5,2)=0.5D+00*(Y(4,2)+Y(6,2))
      Y(MAXXMM-1,2)=Y(5,2)+PITCH
      Z(2,2)=2.0D+00*Z(2,3)-Z(2,4)
      Z(MAXX,2)=Z(2,2)
      Z(3,2)=Z(2,2)
      Z(MAXXM,2)=Z(2,2)
      Z(4,2)=Z(2,2)
      Z(MAXXMM,2)=Z(2,2)
      Z(5,2)=0.5D+00*(Z(4,2)+Z(6,2))
      Z(MAXXMM-1,2)=Z(5,2)
C------ VALUES OF COORDINATES IN IMAGINARY ROWS
      DO 60 I=1,IMID
      II=MAXX2-I
      X(I,1)=X(II,3)
      X(II,1)=X(I,3)
      Y(I,1)=Y(II,3)-PITCH
      Y(II,1)=Y(I,3)+PITCH
      Z(I,1)=Z(II,3)
      Z(II,1)=Z(I,3)
60    CONTINUE
C------ VALUES OF COORDINATES IN IMAGINARY COLUMNS
      DO 70 J=1,MAXY
      X(2,J)=0.5D+00*(X(2,J)+X(MAXX,J))
      Y(2,J)=0.5D+00*(Y(2,J)+Y(MAXX,J))
      Z(2,J)=0.5D+00*(Z(2,J)+Z(MAXX,J))
      X(MAXX,J)=X(2,J)
      Y(MAXX,J)=Y(2,J)
      Z(MAXX,J)=Z(2,J)
      X(1,J)=X(MAXXM,J)
      Y(1,J)=Y(MAXXM,J)
      Z(1,J)=Z(MAXXM,J)
      X(MAXXP,J)=X(3,J)
      Y(MAXXP,J)=Y(3,J)
      Z(MAXXP,J)=Z(3,J)
70    CONTINUE
C------ IMAGINARY POINTS AT UPSTREAM INFINITY
      X(IMID,1)=2.*X(IMID,2)-X(IMID,3)
      X(ILHS,1)=X(IMID,1)
      X(IRHS,1)=X(IMID,1)
      Y(IMID,1)=Y(IMID,2)
      Y(ILHS,1)=Y(ILHS,2)
      Y(IRHS,1)=Y(IRHS,2)
      Z(IMID,1)=Z(IMID,2)
      Z(ILHS,1)=Z(IMID,2)
      Z(IRHS,1)=Z(IMID,2)
C------ IMAGINARY POINTS AT DOWNSTREAM INFINITY
      X(2,1)=2.*X(2,2)-X(2,3)
      Y(2,1)=Y(2,2)
      Z(2,1)=Z(2,2)
      X(MAXX,1)=X(2,1)
      Y(MAXX,1)=Y(2,1)
      Z(MAXX,1)=Z(2,1)
      X(MAXXM,1)=X(2,1)
      Y(MAXXM,1)=Y(MAXXM,2)
      Z(MAXXM,1)=Z(MAXXM,2)
      X(3,1)=X(2,1)
      Y(3,1)=Y(3,2)
      Z(3,1)=Z(3,2)
      X(1,1)=X(MAXXM,1)
      Y(1,1)=Y(MAXXM,1)
      Z(1,1)=Z(MAXXM,1)
      X(MAXXP,1)=X(3,1)
      Y(MAXXP,1)=Y(3,1)
      Z(MAXXP,1)=Z(3,1)
      RETURN
      END
      SUBROUTINE SURF(K,CALE)
C************************************************************
C*    ROUTINE WHICH DETERMINES THE INTERSECTION CURVES      *
C*      BETWEEN THE BLADE AND THE AXISYMM. SURFACES         *
C************************************************************
      COMMON/BLK1/ HR(200),VR(200),XHUB(51),RHUB(51),XDUCT(51)
     1            ,RDUCT(51),XCELL,YCELL,ZCELL,PI,PITCH,MAXP,IMID,ILHS
     2            ,IRHS,MAXXM,MAXX,MAXXP,MAXYM,MAXY,MAXYP,MAXZM,MAXZ
     3            ,MAXZP,KTIP,MHUB,MDUCT
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      COMMON/BLK3/ X(200,51),Y(200,51),Z(200,51)
      COMMON/BLK4/ XP(200,51),YP(200,51),ZP(200,51)
      DOUBLE PRECISION HR,VR,XHUB,RHUB,XDUCT,RDUCT,PI,PITCH,S,F,SI,FI,
     1          RAD,RH,XCELL,YCELL,ZCELL,XP,YP,R,RR,FP,FPP,FPPP,CALE
      IMUN=MAXX+2
      IMAM=MAXX-2
      R=RAD(K)
C------ RHUB AND RDUCT SHOULD BE EVALUATED AT THE SAME X-PNT
      DO 100 I=1,MDUCT
      S(I)=XDUCT(I)
      F(I)=RDUCT(I)
100   CONTINUE
      DO 110 I=1,MHUB
      SI(I)=XHUB(I)
110   CONTINUE
      CALL SPLIF(1,MDUCT)
      CALL INTPL(1,MHUB,1,MDUCT)
      DO 120 I=1,MHUB
      HR(I)=FI(I)
120   CONTINUE
C------ SHAPE OF EACH COAXIAL AXISYMMETRIC R=R(X) SURFACE
      RR=R-CALE
      DO 130 I=1,MHUB
      S(I)=XHUB(I)
      F(I)=RHUB(I)+RR*(HR(I)-RHUB(I))
130   CONTINUE
      CALL SPLIF(1,MHUB)
C------ THE LARGE LOOP STARTS HERE
      DO 1000 J=2,MAXY
C------ DETERMINATION OF THE MOST LEFTWARD AND RIGHTWARD
C------ POINT ON THE EACH ELLIPSE J=CONST.
      IF(J .EQ. 2) GO TO 800
      MXM=IMID-6
      MXP=IMID+6
      DO 400 I=MXM,MXP
      IM=I-1
      IF(X(I,J) .GE. X(IM,J)) GO TO 405
400   CONTINUE
405   CONTINUE
      ILEFT=I-1
      IRIGHT=MAXX
      MXM=MAXX-5
      DO 410 I=MXM,MAXX
      II=IMUN-I
      IIP=II+1
      IM=I-1
      IF(X(I,J) .LT. X(IM,J)) GO TO 415
      IF(X(II,J) .LT. X(IIP,J)) GO TO 420
410   CONTINUE
      GO TO 425
415   CONTINUE
      IRIGHT=I-1
      GO TO 425
420   CONTINUE
      IRIGHT=IIP
425   CONTINUE
      IF(IRIGHT .EQ. 2) IRIGHT=MAXX
      IF(IRIGHT .EQ. MAXX) GO TO 600
      IF(IRIGHT .GT. 10) GO TO 500
      DO 450 I=IRIGHT,ILEFT
      SI(I)=X(I,J)
450   CONTINUE
      CALL INTPL(IRIGHT,ILEFT,1,MHUB)
      DO 455 I=IRIGHT,ILEFT
      Z(I,J)=FI(I)
455   CONTINUE
      DO 460 I=ILEFT,MAXX
      SI(I)=X(I,J)
460   CONTINUE
      DO 465 I=2,IRIGHT
      II=IMAM+I
      SI(II)=X(I,J)
465   CONTINUE
      IRR=IMAM+IRIGHT
      CALL INTPL(ILEFT,IRR,1,MHUB)
      DO 470 I=ILEFT,MAXX
      Z(I,J)=FI(I)
470   CONTINUE
      DO 480 I=2,IRIGHT
      II=IMAM+I
      Z(I,J)=FI(II)
480   CONTINUE
      GO TO 900
C------ CASE WHEN (IRIGHT .GT. 10)
500   CONTINUE
      DO 510 I=ILEFT,IRIGHT
      SI(I)=X(I,J)
510   CONTINUE
      CALL INTPL(ILEFT,IRIGHT,1,MHUB)
      DO 520 I=ILEFT,IRIGHT
      Z(I,J)=FI(I)
520   CONTINUE
      DO 530 I=IRIGHT,MAXX
      SI(I)=X(I,J)
530   CONTINUE
      DO 540 I=2,ILEFT
      II=IMAM+I
      SI(II)=X(I,J)
540   CONTINUE
      ILL=IMAM+ILEFT
      CALL INTPL(2,ILL,1,MHUB)
      DO 550 I=IRIGHT,MAXX
      Z(I,J)=FI(I)
550   CONTINUE
      DO 560 I=2,ILEFT
      II=IMAM+I
      Z(I,J)=FI(II)
560   CONTINUE
      GO TO 900
C------ CASE WHEN (IRIGHT .EQ. MAXX)
600   CONTINUE
      DO 650 I=2,MAXX
      SI(I)=X(I,J)
650   CONTINUE
      CALL INTPL(2,ILEFT,1,MHUB)
      CALL INTPL(ILEFT,MAXX,1,MHUB)
      DO 700 I=2,MAXX
      Z(I,J)=FI(I)
700   CONTINUE
      GO TO 900
800   CONTINUE
      X(IMID,2)=X(ILHS,2)-1.
      DO 825 I=3,MAXXM
      SI(I)=X(I,2)
825   CONTINUE
      CALL INTPL(3,ILHS,1,MHUB)
      CALL INTPL(IRHS,MAXXM,1,MHUB)
      DO 850 I=3,MAXXM
      Z(I,2)=FI(I)
850   CONTINUE
900   CONTINUE
1000  CONTINUE
      RETURN
      END
      SUBROUTINE PISMO(K,N1,LXYZ)
C**************************************************************
C*ROUTINE WHICH WRITES X-Y-Z COORDINATES (LXYZ=1) OR X-THETA-R*
C*   COORDINATES (LXYZ=0) IN UNFORMATED FORM ON UNIT NO.N1    *
C**************************************************************
      COMMON/BLK1/ HR(200),VR(200),XHUB(51),RHUB(51),XDUCT(51)
     1            ,RDUCT(51),XCELL,YCELL,ZCELL,PI,PITCH,MAXP,IMID,ILHS
     2            ,IRHS,MAXXM,MAXX,MAXXP,MAXYM,MAXY,MAXYP,MAXZM,MAXZ
     3            ,MAXZP,KTIP,MHUB,MDUCT
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      COMMON/BLK3/ X(200,51),Y(200,51),Z(200,51)
      COMMON/BLK4/ XP(200,51),YP(200,51),ZP(200,51)
      DIMENSION YCAR(200,51)
      DOUBLE PRECISION HR,VR,XHUB,RHUB,XDUCT,RDUCT,PI,PITCH,XCELL,TH,
     1XP,YP,YCELL,ZCELL,DSIN,DCOS,S,F,SI,FI,FP,FPP,FPPP,RAD,RH,RO
      MAXX2=MAXXP+1
      IF(LXYZ .EQ. 0) GO TO 185
C------ (X,Y,Z) COORDINATES FOR OUTPUT
      DO 180 J=1,MAXY
      DO 180 I=1,MAXXP
      TH=Y(I,J)
      RO=Z(I,J)
      YCAR(I,J)=RO*DSIN(TH)
      ZP(I,J)=RO*DCOS(TH)
180   CONTINUE
      GO TO 195
185   CONTINUE
C------ (X,THETA,R) COORDINATES FOR OUTPUT
      DO 190 J=1,MAXY
      DO 190 I=1,MAXXP
      YCAR(I,J)=Y(I,J)
      ZP(I,J)=Z(I,J)
190   CONTINUE
195   CONTINUE
      DO 240 I=1,MAXXP
      X(I,MAXYP)=2.0D+00*X(I,MAXY)-X(I,MAXYM)
      YCAR(I,MAXYP)=2.0D+00*YCAR(I,MAXY)-YCAR(I,MAXYM)
      ZP(I,MAXYP)=2.0D+00*ZP(I,MAXY)-ZP(I,MAXYM)
240   CONTINUE
      DO 250 I=1,IMID
      II=MAXX2-I
      X(I,2)=0.5D+00*(X(I,2)+X(II,2))
      X(II,2)=X(I,2)
250   CONTINUE
      WRITE(N1) ((X(I,J),YCAR(I,J),ZP(I,J),I=1,MAXXP),J=1,MAXYP)
      RETURN
      END
      SUBROUTINE SPLIF(M,N)
C*************************************************************
C*     MODIFIED ROUTINE FROM: SUPERCRITICAL WING SECTIONS II *
C*     CUBIC SPLINE IN ONE SPATIAL COORD. IS FITTED TO DATA  *
C*     ARRAY "F" AT NODES "S" FROM INDEX "M" TO INDEX"N"     *
C*************************************************************
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      DOUBLE PRECISION S,F,SI,FI,FP,FPP,FPPP,FPPPP,B,C,D,DS,DF,FF
     1 ,RAD,RH,U,V
      K=IABS(N-M)
      K=(N-M)/K
      I=M
      J=M+K
      DS=S(J)-S(I)
      D=DS
      IF (DS) 11,81,11
11    DF=(F(J)-F(I))/DS
      U=-1.0D+00
      V=0.0D+00
      GO TO 25
21    I=J
      J=J+K
      DS=S(J)-S(I)
      IF (D*DS) 81,81,23
23    DF=(F(J)-F(I))/DS
      B=1.0D+00/(DS+DS+U)
      U=B*DS
      V=B*(6.0D+00*DF-V)
25    FP(I)=U
      FPP(I)=V
      U=(2.0D+00-U)*DS
      V=6.0D+00*DF+DS*V
      IF((J-N) .NE. 0) GO TO 21
      V=FPP(I)/(1.0D+00+FP(I))
      B=V
41    DS=S(J)-S(I)
      U=FPP(I)-FP(I)*V
      FPPP(I)=(V-U)/DS
      FPP(I)=U
      FP(I)=(F(J)-F(I))/DS-DS*(V+U+U)/6.0D+00
      V=U
      J=I
      I=I-K
      IF(J-M) 41,51,41
51    I=N-K
      FPPP(N)=FPPP(I)
      FPP(N)=B
      FP(N)=DF+D*(FPP(I)+B+B)/6.0D+00
81    RETURN
      END
      SUBROUTINE INTPL(MI,NI,M,N)
C*************************************************************
C*     MODIFIED ROUTINE FROM: SUPERCRITICAL WING SECTIONS II *
C*         INTERPOLATION USING PIECEWISE TAYLOR SERIES       *
C*************************************************************
      COMMON/BLK2/ S(200),F(200),SI(200),FI(200),FP(200),FPP(200),
     1             FPPP(200),RAD(200),RH(200),RAD1(51),RAD2(51)
      DOUBLE PRECISION S,F,SI,FI,FP,FPP,FPPP,FPPPP,C,SS,FF,D,RAD,RH
      K=IABS(N-M)
      K=(N-M)/K
      I=M
      MIN=MI
      NIN=NI
      D=S(N)-S(M)
      IF (D*(SI(NI)-SI(MI))) 11,13,13
11    MIN=NI
      NIN=MI
13    KI=IABS(NIN-MIN)
      IF (KI) 21,21,15
15    KI= (NIN-MIN)/KI
21    II=MIN-KI
      C=0.0D+00
31    II=II+KI
      SS=SI(II)*1.0D+00
33    I=I+K
      IF (I-N) 35,37,35
35    IF (D*(S(I)-SS)) 33,33,37
37    J=I
      I=I-K
      SS=SS-S(I)
      FPPPP=C*(FPPP(J)-FPPP(I))/(S(J)-S(I))
      FF=FPPP(I)+0.25D+00*SS*FPPPP
      FF=FPP(I)+SS*FF/3.0D+00
      FF=FP(I)+0.5D+00*SS*FF
      FI(II)=F(I)+SS*FF
      IF (II-NIN) 31,41,31
41    RETURN
      END
